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Abstract. We study the convergence of multigrid schemes for the Helmholtz equation, 
focusing in particular on the choice of the coarse scale operators. Assuming that the 
coarse scale solutions should approximate the true solutions, the coarsest level grid must 
have a minimum number, say G c , of points per wavelength. The oscillating nature of the 
solutions implies the requirement G c > 2. However, in examples the requirement can be 
more like G c > 8, in a trade-off involving also the amount of damping present and the 

>— ^ number of multigrid iterations. We conjecture that the main factor affecting G c is the 

difference in phase speeds between the coarse and fine scale operators. Standard 5-point 
finite differences in 2-D are our first example. A new coarse scale 9-point operator is 
constructed to match the fine scale phase speeds. We then compare phase speeds and 
multigrid performance of standard schemes, with a new scheme that uses the new coarse 
scale operators and a modified smoother. The required G c is reduced from about 8 to 
about 4, while at the same time less damping is needed so that waves propagate over ~ 
fOO wavelengths in the new scheme. Next we consider extensions of the method to more 
general cases. In 3-D comparable results are obtained with standard 7-point differences 

j / - , and optimized f 5-point coarse grid operators, leading to an order of magnitude reduction 

in the number of unknowns for the coarsest scale linear system. Finally we show how 
to include PML boundary layers, using a regular grid finite element method. Matching 
coare scale operators can easily be constructed for other discretizations. The method 
is therefore potentially useful for a large class of discretized high-frequency Helmholtz 

i— ^ equations. 
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1. Introduction 



The high-frequency Helmholtz equation forms a major challenge in numerical analysis. 

■^J- When the domain is large compared to the wavelength, i.e. large k, the large indefinite 

■^ systems resulting from the discretization are difficult to solve. High-frequency Helmholtz 

problems have applications in various simulation problems and in inverse problems, e.g. 

in exploration seismology and acoustic scattering. 

Multigrid methods have been used for Helmholtz equations, but generally not in the 
conventional setup as used in elliptic problems (see e.g. |16|). The reason is that con- 
ventional multigrid methods perform poorly Issues are the divergence of the iterative 
smoothing methods and the resolution of the coarse grid. To be precise, denote by G c the 
minimum number of grid points per wavelength in the coarsest grid. If the coarse problem 
solutions are to approximate the solutions of the original problem, the oscillatory nature 
of the solutions leads to the requirement G c > 2. However, in examples the requirement 
is more like G c > 8. 

For large k problems the coarsest scale problem can still be quite large, leading to large 
memory requirements. Shifted Laplacian methods [7] provide a way to address this. They 
use a multigrid cycle as a preconditioner for an outer iterative method like BiCGSTAB. 
However, the multigrid is not applied to the original operator, but to a modified operator 
with a strong damping term added. Coarse grids coarser than the wavelength scale are 
used (G c < 1) (see also [5]). This leads to a method that is convergent and with relatively 
low memory use and cost per iteration, but with quite large iteration numbers, increasing 
with the domain size. 

i 
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Direct methods on the other hand also remain relevant, in particular in view of recent 
results on double sweeping domain decomposition methods, like the moving PML method 
of [6] or the double sweeping method of [H]. The double sweeping domain decomposition 
methods appear to have the best, near- linear scaling for the cost per solve, but have the 
disadvantage of a large memory use. See [12J for a parallel implementation of the method 
of [6] . Other methods include those in J2J [3], |T7] . In all cases it is important to distinguish 
the behavior in non-resonant cases vs. resonant cases. In variable coefficient media with 
resonances the performance of the iterative methods discussed above tends to deteriorate 
strongly. 

In this paper we study multigrid methods from the following point of view. We remain 
in the regime G c > 2, aiming at fast convergence, say < 20 iterations for reduction of the 
residual by 10 -6 . However, we challenge the conventional choices made in the design of 
multigrid methods. With new methods we obtain a reduction of G c from about 8 to about 
4. As examples we study the standard 5-point and 7-point stencils (in 2-D and 3-D resp.) 
for the original, fine scale operator, and a regular grid finite element method for which we 
show how to include PML boundary layers. 

To obtain these results we briefly analyze and modify the smoothing steps in the al- 
gorithm, and then focus on the choice of the coarse grid operators. Standard choices for 
coarse scale operators are the Galerkin approximation, or to use the same discretization 
scheme as the fine scale operator. We will show that these are suboptimal by constructing 
optimized coarse scale operators that perform much better, yielding the aforementioned 
factor 2 improvement of G c in combination with substantially less damping required in the 
equation. Note that a factor 2 improvement in G c reduces the number of unknowns for the 
coarsest scale direct solver by 8 in 3-D, and the computational cost and memory use by 
an even higher factor. Moreover, there is nothing special about the discretizations studied 
in this paper, and we conjecture that in fact our method is applicable to any regular grid 
fine scale discretization. 

We now introduce the setup in more detail. The Helmholtz equation reads 

(1) Lu d = -Au-((l + ai)k{x)) 2 u = f, inficT 

where A is the Laplacian and is a rectangular block. We assume Dirichlet boundary 
conditions the boundary d£l. PML layers will be present in some of the examples, they 
will be discussed in section [5T2| The parameter k in general depends on i £ U. Here a is 
a parameter for the damping, that will mostly be independent of qj The corresponding 
undamped operator will be denoted by 

(2) H = -A-k 2 . 

The imaginary contribution iak leads to exponential decaying solutions. For example, in 
1-D for constant k there are solutions 

fo\ ikx—akx 

This limits the propagation distance of the waves. This distance, measured by the number 
of wavelengths for the amplitude to be reduced by a factor 10, is given by 

Values of a will range from 1.25 -10 3 to 0.02, small enough for applications. In presence 
of PML layers a will be set to zero. The experiments will show that, for the optimized 
coarse grid methods good convergence starts roughly around 0.0025 or at D(a) ~ 145 
wave lengths assuming around 4 points per wavelength in the coarsest grid. 



Other authors sometimes use —A — (1 + ia)k 2 , i.e. with the factor (1 + ia) outside the square. The 
sign of ia is related our choice of temporal Fourier transform /(£) = e _I "'/(o;). 



MULTIGRID FOR THE HELMHOLTZ OPERATOR 3 

Multigrid methods consist of several components. For example, a two-grid method con- 
sists of the following steps: presmoothing using an iterative smoothing method; restriction 
of the residual to the coarse grid; solving a coarse grid equation; prolongation of the so- 
lution to the fine grid; suppressing remaining errors using postsmoothing, using again an 
iterative smoothing method. In our method we introduce two main changes compared to 
conventional multigrid methods. The first concerns the smoothing, the second the coarse 
grid operator. In a two-grid method, this will lead to three discrete versions of the oper- 
ator L: the original, fine scale operator Lh, the operator used for smoothing L$ t h and a 
coarse scale operator L c ,2h- 

As for the smoothing, it is well known that iterative smoothing methods are unstable 
for Helmholtz problems. We find that this can be solved by using instead the Laplacian 
operator in the smoother, or, when a PML layer is present, by using the positive definite 
Helmholtz operator —A + k 2 . The error caused by this modification is handled by the 
outer iterative method, in such a way that the overall effect of this modification is positive. 

The second component to be modified is the coarse grid operator. Our claim is that the 
coarse problem should have the same phase speed or numerical dispersion as the fine scale 
problem. Indeed over large distances the phase errors lead to large differences between 
the approximate, coarse grid solution and the true solution. Alternatively, from Fourier 
analysis of multigrid methods one can argue that the inverse of the symbols for the coarse 
scale and original operators should match, which in turn implies that the zero-sets and 
hence the phase speeds should match. To obtain coarse scale operators with these proper- 
ties we optimize the coefficients in the finite difference scheme, minimizing the difference 
of the symbols, with a strong weighting around the zeroset of the fine scale symbol. 

We first take the example of standard 5-point finite differences in 2-D. After computing 
the optimized coefficients, we compare the phase speeds associated with various finite 
difference discretizations, analyze the two- grid convergence factors for the conventional 
and modified methods, and perform numerical experiments. The multigrid convergence is 
good or poor, precisely when the phase speeds match well or poorly, respectively. Based 
on the two-grid convergence factors, this conclusion is derived in constant coefficient, 2-D 
case. Numerical experiments show that similar conclusions hold for the multigrid method, 
when the coefficient k is spatially varying. In addition we study in 2-D the optimized finite 
difference scheme by Jo, Shin and Suh [9]. Their scheme has accurate phase speeds for all 
values of G (the number of grid points per wavelength) for G down to 4, which means the 
phase speeds match well if this scheme is used at all levels of the multigrid method. It 
also performs well. The advantage of our own optimized scheme is that it can be adapted 
for any regular grid discretization of the original problem. 

Quantitative results are as follows. We observe that for a in the order of 0.25 % we still 
have good convergence for 4 points per wave length with the optimized method. Regular 
multigrid with 8 points per wave length however, only performs well if a > 2 %. Note that 
the actual value of G c , whether it is say 8 or more, or around 3.5 or 4, can make a great 
difference. Reducing the grid by a factor of 2 in each direction, leads in 3-D to a 2 3 fold 
reduction in the number of unknowns and an even larger reduction in the computational 
cost of the sparse factorization and the memory use, up to 2 6 for the sparse factorization 
and 2 9 ' 2 for the memory use [8]. 

Next we consider more general settings than 5-point finite differences. For the 3-D 
problem we construct optimized coarse scale operators for the standard 7-point finite 
difference scheme. Numerical tests show similar results as for the 2-D case, with good 
convergence for G c > 4. 

In practice problems on rectangular domains often occur in combination with absorbing 
boundary conditions or absorbing boundary layers, such as PML layers [H Uj. These 
then provide damping in the equation, so that we set a = 0. In numerical experiments, 
straightforward insertion of a PML layer was observed to lead to very poor convergence 
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or no convergence at all in the situations we are interested in (G c > 4). To address this 
we consider an adapted coarsening strategy. In the PML layer no coarsening takes place 
in the direction normal to the boundary. This is the direction in which a fast variation 
of the coefficient a used in the PML takes place. Because this leads to (partly) irregular 
grids, it is natural to consider this in the context of finite elements. Optimized coarse grid 
operators for a first order rectangular finite element discretization are computed, and we 
show in 2-D examples that this again results in good convergence of the multigrid method 
for G c > 4, and propapagation distances of up to 200 wavelengths. 

The setup of the paper is as follows. The next section focusses on the phase speeds. 
Here we construct the new, optimized finite difference operators for the various cases 
and compare the phase speed errors in the standard and the new coarse grid operators. 
Then in section [3] we describe the multigrid methods used in this study. In section [4] 
we present Fourier analysis of the two- grid methods. Section [5] contains the results of 
numerical simulations: First the comparison of the standard and the new methods for 2-D 
finite differences, and then the extensions to multigrid, to 3-D and to problems with PML 
boundary layers. 

2. Phase speeds and optimized coarse scale operators 

We first recall the notions of dispersion relation and phase speed for constant coefficient 
linear time dependent partial differential equations |15| . If P is such an operator, and 

(5) P(f , u) = e -^+^(iV^-^*), 

denotes its symbol, then the dispersion relation is the set of (£, to) where P(£, oS) = 0. 
For the Helmholtz equation, the symbol is a function of the spatial wave number £, with 
parameter w 

(6) H^-uj) = e-^ x {He^ x ) 

and the dispersion relation is understood to be the w-dependent set 

(7) {£€S|ff(£w)=0}, 

where for the continuous operator S = W 1 , and it is assumed that a = 0. For £ such that 
H(£; uj) = 0, the number t£L is called the phase speed v p ^ associated with a plane wave 
solution. For the continuous operator v p h is constant equal to c. 

For finite difference discretizations of the Helmholtz operator the definitions (pi) and 
(I7| remain valid except that S is given by the fundamental domain S = [—Tr/h,7r/h] n . In 
the typical case that H is increasing along half lines from the origin, the phase speed is 
a function of angle. We will compute the dimensionless phase speed -^ as a function of 
another dimensionless quantity, the number of points per wavelength G, or its inverse 1/G 
(G is related to the dimensionless quantity kh by kh = ^j). 

This discussion will mostly involve only two scales, a fine scale and a coarse scale with 
double the grid parameter. In multigrid with more than two levels, the additional levels 
are assumed to be increasingly fine, because the coarsest level is restricted by the wave 
length. In general we expect that the phase speed difference between the two coarsest 
levels is the most important, since for finer grids these differences automatically become 
smaller as the discretization becomes a better approximation of the true operator. 

Next we first show the behavior of the phase speeds for some well known finite differ- 
ence schemes and then construct a finite difference method that is optimized so that its 
dispersion relation matches that of the standard 5-point method. After this we consider 
the 7-point operator in 3-D and the finite element method used with PML layers in 2-D. 
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2.1. Phase speeds of some well known finite difference schemes. The standard 
5-point finite difference discretization (fd5) of the 2-D Helmholtz operator is given by 

(8) (H^ 5 u)i tj = h~ 2 (4uij - Ui-i,j - iM+ij - iHj-i - u iJ+1 ) - k 2 mj. 
Its symbol is easily shown to be 

(9) ajp(£) = /r 2 (4 - 2cos(/tCi) - 2cos(h&)) - k 2 , 

where £ = (£i, £2) denotes the wave vector. The phase speed as a function of angle is easily 
computed using a root finding algorithm and shown in Figure [TVa) . The phase speeds as 
a function of |£| are given for 4 angles 0°, 15°, 30° and 45°. 

When such a scheme is used in a multigrid method, say a two-grid method for the 
purpose of this argument, standard choices for the coarse scale operator are to use the 
same discretization, or to use a Galerkin discretization. The comparison between the phase 
speeds of a coarse scale operator fl^ and a fine scale operator H^?2 ^ s gi ven m Figure l[b). 
Here G refers to the coarse scale. The Galerkin method, using "full weighting" restriction 
and prolongation operators |16| is easily shown to have symbol 



of® = 3h- 2 - -k 2 + (-h- 2 - -fe 2 )(cos(£i/0 + cosfe/i)) 
(10) 16 16 

+ (" 2 ^ 2 - - A; 2 )(cos((6 + &)h) + cos((6 - &)fc)). 

The phase speeds and the phase speed difference between a coarse scale Galerkin and a 
fine scale fd5 method are shown in Figure [TVc) and (d). The phase speed differences are 
on the order of 0.01 or 0.02 for G = 8 and larger for G smaller. 

We contrast this with the optimized finite difference method of Jo, Shin and Suh [9j, 
sometimes called the mixed grid operator. This operator, here called JSS operator for 
brevity, is given by 

(11) 

(H^ s u)ij = ((2 + 2a)bT 2 - ck 2 )u iy j + (-ah~ 2 - dk 2 )(u i+ltj + m-ij + Ujj + i + Uij-t) 

1 1 — a r-2 (1 — c — 4cf) 2 
+ ( ^—n r k )(Uj + ij + i + Uj-ij+i + u»+ij_i + Ui_ij_i) 

where a = 0.5461, c = 0.6248 and d = 0.9381 • 10 _1 . The symbol of this operator hence 
equals 

ajs s (£) = (h- 2 (2 + 2a) - k 2 c) + 2(-h~ 2 a - ^(cos^i/i) + cosfah)) 

+ 2{-h- 21 -^- - (1 ~ C 4 ~ 4d) fe 2 )(cos((ei + &)h) + cos((6 - £ 2 )h)). 

The phase speed and relative phase speed difference between H 3 ^ and H 3 ^ s are given in 
Figure fife) and (f). The phase speed difference between coarse and fine level operators, 
both using the JSS method, is reduced to around 2.5 • 10~ 3 for G down to 4, i.e. a very 
substantial improvement both in G and in the size of the phase speed difference. 

The operators L^ 5 , L| a and L J ^ S , with a/0 are obtained simply by replacing k 2 by 
((l-Ha)fc) 2 . 

2.2. An optimized finite difference scheme. We now construct a finite difference 
method with phase speed matching that of the standard 5-point method, with compara- 
ble accuracy as the JSS method, and even down to G = 3. The construction is easily 
generalized to 3-D and to other fine scale methods. 
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Figure 1. Dimensionless phase speed curves and difference between fine 
and coarse scale phase speeds; (a), (c) and (e): phase speed for the fd5, 
Galerkin and JSS method; (b), (d) and (f): phase speed differences vf^ 5 — 

& v f l ~ "81 > and v h SS - v h/2> with an S les 0°, 15°, 30°, 45° 



J h/2i 



We consider finite difference operators of the following form 

H^ p u(xij) = - k 2 u itj + h~ 2 (a\Uij + a 2 (uj_ij + u i+ ij + Uij-i + Ujj+i) 
(13) + a 3 (n i _i J -_i + Ui_ij + i + u i+ ij-i + u i+ i : j+i) 

+ a 4 (Ui-2J + Ui +2 ,j + Ui,j-2 + Ui,j+2)), 

where we include the possibility that a^ = 0, or a^ = 0. Partly based on the numerical 
results below, we will in fact take 04 = in most of work. The associated symbol is 



(14) 



a 



opt 



- k 2 + h (ai + 2a 2 (cos(/i£i) + cos(/i£ 2 )) 

+ 2a 3 (cos(/i(£i + 6)) + cos(/i(fi - &)) + 2a 4 (cos(2/i£i) + cos(2/i&)). 



-opt 



For the operator L fe p with a/0, the a,j are still a function of k, and the term k' 
replaced by ((1 + ai)k) 2 Uij. 



UiJ is 
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The coefficients a,j, will be chosen as a function of (/i, /if, k), where h = h c is the grid 
distance of the grid the operator is used in, i.e. the coarse grid, and /if represents the 
grid distance of the fine scale operator to be approximated. However, from dimensional 
considerations, the coefficients Oj only have to depend on two parameters, namely aj = 
aj(-^,hk), where the parameter hk can also be replaced by the parameter 1/G as noted 
above. 

The optimization procedure implements two main features 

• The optimized coarse symbol is an approximation to the fine symbol for all £. 

• However, the accuracy of the approximation depends on the position in the £ plane. 
In particular, the approximation should be highly accurate around the zero-set of 
the fine symbol in order to reproduce the phase speeds. 

This is done in an ad hoc fashion, choosing a set of points Pj in the £ plane, with for each 
point an associated weights Wj > and subsequently finding (01, 02, ^3, 04), (or (01, 02, 03) 
in case 04 = or (ai, 02, 04) in case 03 = 0) that minimize 



(15) 



E^Ca^i)-^^)) 5 



The starting point for the choice of the Pj is the closed curve r 
nates), such that 

' fine (r o (0)cos0,r o (0)sin0) 



ro(9) (in polar coordi- 



a ' 



0. 



(16) 

We choose ng points on the set r = r(9), equally spaced in angle, with Wj = 1. Then we 
choose points at 0.85r(6>) and at 1.15r(0), with weight 1/6, and points at O.5r(0) and at 
1.6r(0) with weight 0.05. Finally we add points at (0,0), (0, 7r), (it, it) with weights 0.01, 
0.01, 0.01. These ad hoc choices are justified by the good behavior of the obtained phase 
speeds, as well as the obtained multigrid methods (see below). But other choices could be 
possible as well. 

We computed results for 3 coefficients, 0,1,0,2,(13. For hf/h = 0.5 the result is given in 
Figure [2] We observe that the numerical optimization indeed produces a result and the 
coefficients appear to depend smoothly on hk. 

Choosing values of hf/h different from 1/2 is useful in multigrid with more than two 
levels, see below. The maximum over 9 and the root-mean-square over 9 of the phase error 
are given in Figure El where hf/h is given by 1/2, 1/4, 1/8 respectively. 

As can be seen from these figures, the phase speed differences are reduced to around 
2.5 • 10~ 3 for G c downto to 3, comparable to JSS, and with in fact a slight improvement 
below G c = 4. 
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Figure 2. Coefficients Oj as function of 1/G = ^; (a) coefficients 01, 02, 03 
when CJ4 = 0; (b) coefficients 01,02,04 when 03 = 0; (c) coefficients 
0,1,02,03,04. 
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Figure 3. (a) Dimensionless phase speed for the opt method; (b) phase 
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Figure 4. Phase speed error for the different values of h{/h versus 1/G. 
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Figure 5. Optimized coefficients and phase speed errors in 3-D; (a), (b) 
the 15-point operator described in the text; (c), (d) the 27-point operator. 
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Figure 6. Optimized coefficients and phase speed errors for 2-D finite elements. 

2.3. Optimized finite difference scheme in 3-D. In 3-D we similarly construct an 
optimized operator to match the standard seven point finite difference discretization of 
the Helmholtz operator. The coefficients of the optimized operator are again denoted 
by 01,02,03,04. These coefficients correspond to contributions from the center, the face- 
centers, the edge-centers and the corners of a 27-point cube, such that the symbol becomes 

T o P t = _ k 2 + h -2( Qi + 2o2 (cos(^i) + cos(o£ 2 ) + cos(/i&)) 

+ 2o 3 (cos(/i(£ 1 + £ 2 )) + cos(o(£i - £2)) + cos(o(£i + £3)) + cos(o(£i - £3)) 

+ cos(/i(£ 2 + &)) + cos(/i(£ 2 - £3))) + 2a 4 (cos(/i(£i + £ 2 + £3)) 

+ cos(/i(fi - £2 + £3)) + cos(o(£i + £ 2 - 6)) + cos(n(6 - 6 - &)))■ 



(7 



(17) 



A very similar optimization method was applied, minimizing a weighted difference of <r opt 
and CT fd7 . Two cases were studied, first the case with 03 = 0, which leads to a 15-point 
operator, and then the case where all four coefficients are allowed to be nonzero. The 
values of the coefficients, and the rms and maximum errors in the phase speed are given 
in Figure [5] 

2.4. Optimized regular grid finite elements. It is straightforward to derive the matrix 
for 2-D regular grid finite elements using first order rectangular (square) elements and 
constant k 2 , see section 5.2 below. One finds 

)i,3 =(o~ q h2k2 ) U i,J + (~o ~ q h2k2 )( U i+lJ + U i-hj + U M + 1 + u i,j-l) 
+ (~ o ~ ^^ 2 ^ 2 )( n i+l,i+l + u i-l,j+l + u i+l,j~l + Ik-lj-l) 



(18) 



(H F h E u n , 



3 36 

Note that H scales differently with h due to the use of the finite element method instead 
of finite differences. It is straightforward to find optimized coarse scale finite difference 
operators in the form 



(19) 



H 



opt.FE 
h 



u(xi t j) = amij + a,2(ui-i,j + Ui+ij + Ujj-i + Ui,j+i) 

+ 03(«i-l,i-l + Wi-lJ+1 + u i+l,j-l + Ui+l,j+l) 



The same optimization scheme is used, except that the points (0,7r), (tt, 7r) are omitted 
from measure of symbol differences. The values of the coefficients, and the max and rms 
errors are given in Figure [6} 



3. The multigrid methods 

The above described operator is used in two-grid and multigrid methods. In this section 
we describe the multigrid methods used. For background on multigrid methods we refer 

to fna. 
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The multigrid method will be used as a preconditioner for a Krylov subspace method 
(GMRES pj). 

We next discuss the choice of smoother. Smoothers apply one or more steps of an 
iterative method for the equation 

(20) L s , h u h = f h . 

This generally involves a choice of the iterative method and of the number of smoothing 
steps (presmoothing and postsmoothing) . Here the smoother will be denoted by S%, v 
being the number of iterations, and will be based on w-Jacobi or successive over relaxation 
(SOR) (which, when the relaxation factor is 1 of course reduces to Gauss Seidel (GS)). 

However, we also consider the operator Ls t h as subject to choice, setting it equal to 
a discretization of —A instead of the Helmholtz operator we are trying to invert. The 
reason to not choose Lg^ equal to the Helmholtz equation has to do with the instability 
in the w-Jacobi and SOR methods for the discrete Helmholtz operator, and will be further 
explained in section [4] on Fourier analysis of the multigrid methods below. Note that this 
type of modification is different from shifted Laplacian methods. In the latter both the 
operator used for smoothing and the operator used for coarse grid correction are modified 
to have additional damping, but they remain equal. 

The choices of grid refinement, and of the restriction and prolongation operators are 
standard. We start with a square grid with grid size h, and apply standard coarsening 
to square grids with grid size 2h, Ah etc. For the restriction operator we mostly use full 
weighting, i.e. it is the tensor product of one-dimensional FW operators of the form, in 
stencil notation, \ [l 2 l] . In 2-D it reads, in stencil notation 

(21) R h = -L 

or, if i, j are coordinates in the coarse grid 

(Rhu)ij = — ( 4u2i,2j + 2(u2i-l,2j + U2i+l,2j + U2i,2j-l + U2i,2j+l) 

(22) V 



"1 


2 


1" 


2 


4 


2 


1 


2 


1 



+ U2i-l,2j-l + U2i+l,2j~l + U2i-l,2j+l + U2i+l,2j+l 

As an alternative, the use of straight injection is briefly studied in section |4j The prolon- 
gation operator is given by 

(23) P h = 4(22,0' • 

The two-grid approximate solution formula is then given by 

(24) M* h = S?K* h S£, 
where 

(25) Kf = I h - P h {L Cy2h )- l R h L h , 

is the coarse grid correction operator. Here L^ and L C 2h denote discretizations of — A — /c 2 , 
that are still to be specified. We focus particularly on the choice of coarse grid operator 
L c ,2h ■ F° r the two-grid method in 2-D the following pairs of fine and coarse grid operators 
will be considered 

(1) Both using standard 5-point finite differences (fd5-fd5) 

(2) Using standard 5-point finite differences at the fine level and a Galerkin approxi- 
mation of this matrix at the coarse level (fd5-gal) 

(3) Using standard 5-point finite differences at the fine level and the new optimized 
operator at the coarse level (fd5-opt) 

(4) Using the operator from [9] at the fine and coarse levels (jss-jss) 
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For multigrid with more than two coarsening levels we consider only the optimized method 
where at each scale the optimized operators is chosen to approximate the fine level finite 
difference operator. This is done using the parameter h Te \ = hfi ne /h collTse , which assumes 
the values 1/2, 1/4, 1/8, etc. In this case the algorithm is based on the V-cycle. 

Comparing the case of two levels with the case of more than two levels, the additional 
levels are added on the fine side, i.e. assuming identical coarsest grid. The coarse level 
behavior is expected to determine the behavior of the algorithm and therefore most em- 
phasis is on the two-grid problem. A few experiments are done using multiple levels, to 
study whether it is indeed the two-grid performance determines the convergence behavior. 

4. Fourier analysis of multigrid methods 

In this section, we compute and analyze two-grid convergence factors by using local 
Fourier analysis as described in [16], chapters 3 and 4. We concentrate on the computation 
of the spectral radius p(M% ) for the two-grid iteration operator M% that determines the 
asymptotic convergence behavior. We consider the 2-D case, using standard 5-point finite 
differences on the fine scale with an optimized finite difference scheme L Q ^h = L^ as well 
as the same five point finite difference operator L c ^h = L^ h for the coarse grid correction 
to show the efficiency of the optimized nine point finite difference scheme L«? . We also 
study the case with the Jo-Shin-Suh operator on the fine and the coarse scale. 

The analysis is done in the dimensionless Fourier domain, i.e. the wave vector is written 
as h~ 1 9 = h~ 1 (9i,$2), so that the dimensionless wave numbers are in [— 7r,7r) 2 (in 2-D). 
Sets of "high" and "low" wave numbers are defined by 

T low = \_n ^\ 2 
(26) L 2' 2/ ' 

T high = [_ 7r>7r )2\ r low_ 

For = (0i, 2 ) € T low we also define 

(O ' O) = (0i,0 2 ), (1,1) = (0i,0 2 ). 

{2,) 0(°' 1 ) = (0 1 ,0- 2 ), 0( 1 '°) = (0- 1 ,0 2 ), 



where 

(28) 



+ vr if 9i < 
- TV if 9i > 0, 



The two-grid operator was given in (24), (25). In local Fourier analysis, for each G x 
we define the four-dimensional space of harmonics by 

(29) E e = span ({ e ^-M | a G {(0, 0), (1, 1), (0, 1), (1, 0)} }) 

The two-grid operator maps this space into itself, and hence can be reduced to a 4 x 4 
matrix for each £ j' low ) which will be denoted by M% h (9). It can be written as 

(30) Ml h {9) = S h {9) U2 Kl h {9)S h {9) u \ 9 e T low , 
where 

(31) Kl h {9) = I h - P h (9)(L c , 2h (29))- 1 R h (9)L h (9), 

where Ih is the 4x4 identity matrix and the meaning of K^ h (9), Ph{9), ^c,2/i(20), Rh(9) 

and Lh(9) will be explained next. The Fourier transformed restriction and prolongation 

operators R^ and P^ are 1x4 and 4x1 matrices respectively, given by 

(32) 

i4 = (i(l + cos6>i)(l + cos0 2 ) ;j(l + cos(9i)(l + cos0 2 ) !(1 + cos0i)(1 + cos0 2 ) |(1 + cos0i)(l + cos0 2 )) 
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and 



(33) P h 



/|(l + cosfli)(l + cos0 2 )\ 

|(l + cos#i)(l + cos6> 2 ) 

|(l + cos#i)(l + cos0 2 ) 
\l{l + cos9 1 ){l + eos0 2 )J 

The operator L c2 /j(2#) is a scalar given by the value of the symbol a(9^ 0,0 ' /h). The 
operator L^ is a diagonal 4x4 matrix given by 

(34) L h {9) = diag(L^(°> )), L h (9^), L h (9^), L h (9^) ), 

in which Lh(9) is the scalar value of the symbol a{9/h). The matrix Sh(9) is also a 4 x 4 
diagonal matrix 

(35) 5 h (0) = diag(^(^°)),5 h (^ 1 ' 1 )),^(^ 1 )),5 h (^ 1 '°))), 

containing the smoothing factors Sh(9), which we now discuss. 

The smoothing factors depend on the iterative method and choice of discretization, see 
[16 1 for their derivation. For standard 5-point finite differences of the Helmholtz operator 
and w-Jacobi they are given by 

(36) S h (0) = 1 - u + , 2 ^ 9L9 (cosfli + cosfl 2 ). 

4 — Arn z 

and for SOR by 

- = (1 - uQ(4 - k 2 h 2 ) + o;(exp(igi) + exp(ifl 2 )) 
{ ' h{ ' 4-A; 2 /i 2 -w(expH(9 1 ) + exp(-i0 2 )) ' 

This concludes the explanation of the computation of M^ h . Note that p\ oc (M^ h ) is 
given by 

(38) Pioc(M 2 h h ) = sup {pioc(Mf (0)) | d G T low } , 

and that p\ oc (M^ h (9)) is the spectral radius of the 4 x 4 matrix Ml h (9). 

4.1. Using the Laplacian for smoothing. Using the smoothing factors just computed, 



the choice of using the Laplacian for smoothing can be motivated. From (36) it follows 
that, for e = (9i,9 2 ) = (0,0) 

(39) 5 fcW = i_ w + _j£_>i 

and hence the iterative smoothing method is unstable in this case. If kh = on the other 
hand the smoothing factor satisfies 

(40) 1 - 2w < S h {9) < 1. 



Similarly observations can be made for ( 37 ) . Hence using a discretization of the operator 
—A in the smoothing steps means the instability disappears (no blow up of errors for 
9 ~ (0,0)). However, this use of a different operator means that the smoothing method 
doesn't converge to the true solution. Instead, in the normalized Fourier domain, it 
converges to 

4-2 cos(fli) - 2 cos(fl 2 ) - k 2 h 2 _ k 2 h 2 

\ / a o — „i a \ o „„„( a \ ~ 



4-2cos(0i)-2cos(# 2 ) 4-2cos(6>i)-2cos(6l 2 ) 

times the solution. The purpose of the smoother is to invert the equation for the high wave 
numbers, as for the low wave numbers the coarse grid correction is used. For (#i,# 2 ) £ 



T^ lg the factor (41) is somewhere between 0.21 and 1 assuming that G c > 4 and hence 



kh < -¥-. Typically this kind of factor is easily dealt with by the outer Krylov method, 
and we conjecture that the stabilizing effect of our choice of Lg /, outweighs the error 
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a 


0.0025 


0.005 


0.01 


(v\,v 2 ) 


(1,1) 


(1,2) 


(2,2) 


(1,1) 


(1,2) 


(2,2) 


(1,1) 


(1,2) 


(2,2) 


SOR1.2 


> 1 


0.4710 


0.2320 


0.7956 


0.2816 


0.2320 


0.5024 


0.2320 


0.2320 


SOR1.1 


> 1 


0.3002 


0.2344 


0.6646 


0.2320 


0.2320 


0.4067 


0.2320 


0.2320 


GS 


> 1 


0.3244 


0.2569 


0.6882 


0.2320 


0.2320 


0.4064 


0.2320 


0.2320 


SOR0.9 


> 1 


0.5014 


0.2871 


0.9235 


0.2966 


0.2320 


0.5356 


0.2330 


0.2320 


SOR0.8 


> 1 


0.9096 


0.3513 


> 1 


0.5249 


0.2327 


0.7678 


0.3083 


0.2320 


SOR0.7 


> 1 


> 1 


0.7272 


> 1 


0.8927 


0.4224 


> 1 


0.5178 


0.2727 


J 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


0.9J 


> 1 


0.5120 


0.4096 


0.8373 


0.5120 


0.4096 


0.6400 


0.5120 


0.4096 


0.8J 


> 1 


0.9487 


0.4054 


> 1 


0.5599 


0.2867 


0.7723 


0.3607 


0.2419 


0.7J 


> 1 


> 1 


0.8041 


> 1 


0.9418 


0.4776 


> 1 


0.5563 


0.3253 


0.6J 


> 1 


> 1 


> 1 


> 1 


> 1 


0.8368 


> 1 


0.8434 


0.4969 


0.5J 


> 1 


> 1 


> 1 


> 1 


> 1 


>1 


> 1 


> 1 


0.7927 



Table 1. Two-grid convergence factors p\ oc (M^ h ), G = 4, fd5-opt. 



introduced. This is confirmed in the numerical simulations below. For the Gauss-Seidel 
method similar behavior occurs. 

4.2. Computation of two-grid convergence factors. The two-grid convergence fac- 



tors are computed using the expression (38), taking the supremum for 6 on regular grid. 



This grid was restricted to the first quadrant because of the symmetry, and taken large 



enough to accurately compute the maximum in ( 38 ) . We have computed the convergence 
factor for three schemes 

(1) The fd5-opt scheme using the five point scheme in the fine grid along with opti- 
mized scheme for the coarse grid 

(2) The fd5-fd5 scheme using the five point scheme on both the fine and the coarse 
grids 

(3) The jss-jss scheme using the JSS operator on both the fine and the coarse grids. 

In Table 111 we present p\ oc (M^ ) and for the fd5-opt scheme, using SOR as smoother 
with u = 0.8, 0.9, 1.0 and 1.1, and with w-Jacobi as a smoother, with u = 0.7, 0.8, 0.9 and 
1.0. In this case the performance for G c = 4 studied. In Table [2] results for the fd5-fd5 
method are given, using the same smoothing methods and G c = 8. Best performance 
occurs roughly for SOR with u = 1 (i.e. the Gauss-Seidel method) and for w-Jacobi with 
uj = 0.8 or 0.9. It is clearly observed that the fd5-opt method requires smaller a for 
convergence, and that it has good convergence factors for G c = 4. For the fd5-fd5 method 
G c = 8 was taken to establish a good choice of smoothing method, as lower values of G c 
gave poor convergence. We have also investigated the SI operator for restriction. However, 
the performance was very poor compared to the full weighting restriction operator. 

In Table [3] we vary both a and the number of grid points per wave length in the coarse 
grid. In this table we also included results for the jss-jss method. We observe that the fd5- 
opt and jss-jss methods are comparable, and much better than fd5-fd5. Concerning the 
number of smoothing steps, for GS (^1,^2) = (1, 1) appears a good choice and (^1,^2) = 
(2, 2) for w-Jacobi. Clearly the damping a plays an important role in the convergence. 
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a 


0.005 


0.01 


0.02 


("1,"2) 


(1.1) 


(1,2) 


(2,2) 


(1,1) 


(1,2) 


(2,2) 


(1,1) 


(1,2) 


(2,2) 


SOR1.1 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.8128 


0.7329 


0.6661 


GS 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.8151 


0.7660 


0.7091 


SOR0.9 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.8364 


0.7932 


0.7469 


SOR0.8 


> 1 


> 1 


> 1 


>1 


> 1 


> 1 


0.8893 


0.8221 


0.7807 


J 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


1.0000 


1.0000 


1.0000 


0.9J 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.9135 


0.8668 


0.8338 


0.8J 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.9568 


0.8864 


0.8498 


0.7 J 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


> 1 


0.9193 


0.8716 



Table 2. 
and G c = 



Two-grid convergence factors p\ oc (M^ h ) for the fd5-fd5 method 





a = 0.00125 


a = 0.005 


a = 0.02 


G = 3.5 


fd5-opt 


fd5-fd5 


jss-jss 


fd5-opt 


fd5-fd5 


jss-jss 


fd5-opt 


fd5-fd5 


jss-jss 


GS 


> 1 




> 1 


0.9147 




> 1 


0.3043 




0.3990 


0.8J 


0.9646 




> 1 


0.3264 




> 1 


0.2901 




0.2986 


0.9J 


0.5722 




> 1 


0.4096 




0.9307 


0.4096 




0.4096 


G = 4 




GS 


> 1 




> 1 


0.6865 




0.7666 


0.2794 




0.3019 


0.8J 


0.4054 




> 1 


0.2866 




0.5265 


0.2320 




0.1572 


0.9J 


0.4096 




> 1 


0.4096 




0.4096 


0.4096 




0.4096 


G = 5 




GS 


> 1 




> 1 


0.4224 




0.5151 


0.2302 




0.2029 


0.8J 


0.6439 




0.8030 


0.2591 




0.2576 


0.1811 




0.1388 


0.9J 


0.6148 




0.7522 


0.4096 




0.4096 


0.4096 




0.4096 


G = 6 




GS 


0.7605 


> 1 


> 1 


0.3044 


> 1 


0.3891 


0.1881 


> 1 


0.1613 


0.8J 


0.7612 


> 1 


0.9803 


0.2577 


> 1 


0.3109 


0.1561 


> 1 


0.1303 


0.9J 


0.7454 


> 1 


> 1 


0.4096 


> 1 


0.4096 


0.4096 


> 1 


0.4096 


G = 7 




GS 


> 1 


> 1 


> 1 


0.3034 


> 1 


0.3446 


0.1797 


1.0058 


0.1613 


0.8J 


> 1 


> 1 


> 1 


0.3131 


> 1 


0.3047 


0.1468 


1.0585 


0.1296 


0.9J 


> 1 


> 1 


> 1 


0.4096 


> 1 


0.4096 


0.4096 


1.0316 


0.4096 


G = 8 




GS 


> 1 


> 1 


> 1 


0.3864 


> 1 


0.2803 


0.1840 


0.8126 


0.1614 


0.8J 


> 1 


> 1 


> 1 


0.3948 


> 1 


0.2819 


0.1472 


0.8456 


0.1296 


0.9J 


> 1 


> 1 


> 1 


0.4096 


> 1 


0.4096 


0.4096 


0.8299 


0.4096 



Table 3. Two-grid convergence factors as a function of 
smoothing method for fd5-opt, fd5-fd5 and jss-jss. 



a. 



G, choice of 



5. Numerical results and discussion 

In this section we compare the convergence of multigrid methods with different choices 
of the coarse scale operators. We start with several experiments for the two-grid method 
in 2-D. Then we investigate the extension to the multigrid method and to the 3-D case, in 
order to establish that the behavior observed in the 2-D, two-grid experiments also occurs 
more generally. 

5.1. The two-grid method in 2-D. For a full comparison we study the trade-off between 
damping present in the equation, the number of grid points per wave length at the coarse 
level G c and the number of iterations required to reduce the residual by a factor of 10 -6 . 
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In 2-D it is relatively easy to perform this kind of experiments, as in 3-D the size of the 
coarse scale problem quickly becomes large for modern desktop machine (a machine with 
8GB memory was used for these experiments using a Matlab implementation). 

We now specify the methods use in these simulations. The regular five point operator 
at the fine scale is combined at the coarse scale with three choices of coarse scale operator, 
namely the regular five point operator, the Galerkin operator and our optimized method. 
These combinations are denoted fd5-fd5, fd5-gal and fd5-opt. We also study the JSS 
operator, used both at the fine and at the coarse scale, referred to as jss-jss. As shown 
below, sharp differences are present between the optimized coarse scale methods on the 
one hand and the fd5 and Galerkin coarse scale methods on the other hand. 

The first series of experiments concerned a constant coefficient medium, i.e. k is con- 
stant equal to (T^, where G c is the number of grid points per wavelength at the coarse 
scale. Experiments were performed with three variants for the smoothing method, first 
the Gauss-Seidel (GS) method, secondly Jacobi relaxation with cj = 0.8 (JO. 8) and thirdly 
the Jacobi relaxation with uj = 0.9. In the first case the number of pre- and postsmoothing 
steps equalled (m, 112) = (1, 1), in the second and third case the choice (m, TI2) = (2, 2) was 
used. We found that all three methods performed comparably for the 5-point and Galerkin 
coarse scale operators. For the JSS method the second and third method performed com- 
parably and better than the first. For the opt method the third method performed mostly 
better than the second, with very few exceptions, and both performed better than GS. 
We therefore display in these experiments only the results for the Jacobi relaxation with 
u = 0.9 and (m,n 2 ) = (2,2). 

The results are given in Table |4j The following can be clearly observed 

• The fd5-fd5 and fd5-gal method require relatively large G c , say > 7 as well as large 
damping of a = 0.01 or 0.02 or larger. 

• The fd5-opt and jss-jss methods perform well with G c as low as 3.5 and a = 0.0025 
or a = 0.005 or larger. 

Here iViter ^ 20 is used as the (somewhat subjective) criterion for good convergence. 
Clearly the optimized methods (jss and opt) perform much better than the conventional 
choices fd5 or gal at the coarse level. 

Next we test this in two variable coefficient media. The first is a random medium, for 
which results are given in Table [5j and the second is the Marmousi model, see Table |6j 
The Marmousi model has size 9200 x 3000 meters, and wave speeds between 1500 and 
5500 ms _1 , see the plot. The value of G c is the minimum value present in the model. 
For the random medium, this minimum is assumed only in small, isolated regions of 
the model, for the Marmousi model it is assumed in the top-layer, i.e. a larger region. 
This perhaps explains the small differences, with the convergence in the random medium 
starting roughly at G c = 3 for the fd5-opt and jss-jss methods. 

5.2. Further experiments: Multigrid, 3-D and PML layers. Further numerical 
experiments are done to establish whether these result extend to multigrid experiments, 
to 3-D and to examples with PML layers. Our method to include PML layers is new and 
is described in this section. 

Results for 2-D multigrid with 2, 3 and 4 levels, using a constant coefficient medium, 
are given in Table [7j In each case the same coarse level grid is used. For small G, the 
convergence deteriorates, but from G c > 4 the multigrid method converges about as fast 
as the two-grid method. For the random medium in Table [5j similar behavior is observed, 
but in this example the multigrid method converges about as fast as the two-grid method 
from G c > 3.5. 

In 3-D the cost and memory use of the sparse factorizations scales worse than in 2-D. 
In this paper we study only a relatively small example of size 80 3 that can be done on a 
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G c 


iVi ter for fd5-fd5,nx 


=511 




iV iter 


for fd5 


-gal,nx= 


=511 




a=1.25e-3 


2.5e-3 0.005 


0.01 


0.02 


a=1.25e-3 


2.5e-3 


0.005 


0.01 


0.02 


6 


MOO 


MOO MOO 


94 


33 


MOO 


MOO 


MOO 


74 


27 


7 


MOO 


MOO MOO 


55 


22 


MOO 


MOO 


MOO 


46 


20 


8 


MOO 


MOO MOO 


36 


16 


MOO 


MOO 


93 


33 


15 


10 


MOO 


MOO 49 


21 


11 


MOO 


MOO 


45 


20 


11 


12 


MOO 


71 28 


14 


9 


MOO 


60 


26 


14 


9 


G c 


iViter 


for fd5-opt,nx= 


=1023 




Mter 


for jss- 


ss,nx= 


1023 




a=1.25e-3 


2.5e-3 0.005 


0.01 


0.02 


a=1.25e-3 


2.5e-3 


0.005 


0.01 


0.02 


3 


MOO 


MOO MOO 


MOO 


MOO 


MOO 


MOO 


MOO 


56 


26 


3.5 


93 


33 19 


15 


13 


MOO 


97 


35 


18 


13 


4 


43 


20 13 


10 


9 


50 


21 


13 


10 


9 


5 


22 


12 9 


7 


7 


26 


13 


9 


7 


7 


6 


23 


12 8 


7 


6 


37 


15 


9 


7 


6 


7 


44 


17 10 


7 


6 


35 


15 


9 


7 


6 


8 


89 


25 12 


8 


6 


26 


13 


8 


6 


5 


10 


94 


25 12 


8 


7 


18 


10 


7 


6 


5 


12 


34 


16 10 


7 


6 


12 


8 


6 


5 


4 



Table 4. Convergence results for the two-grid method in a constant 
medium using different combinations of fine and coarse scale operators. 



regular machine with 8 GB memory using Matlab. This was approximately the largest 
example that could be done in this setup. Two types of optimized operators where used, a 
15-point discretization, using points in the center, on the faces and the corners of a 27-point 
cube, and a 27-point discretization using all points in a 27-point cube. The factorization 
of the 15-point operator used about 30 % less memory than that of the 27-point operator. 
The results on the unit cube with constant wave speed are given in Table [8j We use only 
our optimized method. Generalization of the method of [9j to 3-D exist |11| [3] but we 
have not tested these. 

In 3-D, the observed convergence behavior as a function of a and G c does not differ 
much from the behavior observed in 2-D. 



We next consider the 2-D problem with PML boundary layers. This is important in 
practice - rectangular domains are often encountered in combination with PML boundary 
layers, for example in the seismic problem [11]. As mentioned, we consider a finite element 
discretization for this problem. 

As explained in |10j . in the PML method for the boundary x\ = constant the derivative 
■£- is replaced by , . _ t — -, — sjf—- Here a\ is chosen to be on the internal domain (no 
damping), and increases quadratically from the onset of the damping layer to boundary 
of the computational domain. We introduce PML layers on all four boundaries. It is 
convenient to multiply the equation by (l + iu}~ 1 a\{xij) (l + zuM 02 (X2)), this leads to a 
symmetric operator easily discretized by finite elements. 

Straightforward inclusion of PML layers in the finite difference multigrid method was 
observed to lead to very poor convergence, or no convergence at all. Therefore we consider 
a coarsening strategy where inside the PML layer no coarsening takes place in the direction 
normal to the boundary. The resulting grids are given (schematically) in Figure u\ We 



consider first order rectangular elements on these grids, denoted by 



J 'j- 



The matrix 
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Figure 7. Coarsening strategy to handle PML layers. Inside the PML 
layers (in grey) there is no coarsening in the direction normal to the bound- 
ary. 



elements are given by 



(42) 



Mij-ki 



a\a 2 



dx\ dx\ 






dx 2 dx 2 



k(x) 2 a 1 1 a 2 1 (j) k i(x)(l) k i 



dx, 



where the function <Xj(xj 



contains the modifications for inclusion of the 



1+iOJ 1 (7j(Xj) 

PML layer. 

In the finite element context, the prolongation operator P follows straightforwardly from 
specifying the grids, while the restriction operator R is given by P' . Based on experiments 
with constant coefficient example we choose as a smoother one iteration of the u;-Jacobi 
method for the finite element discretization of — A + k 2 , with relaxation factor uj = 0.8. 

For the coarse grid operators the optimized stencils where computed in section [2] These 
are used in the internal, square grid part of the domain, in such a way that the correspond- 
ing rows of the matrix contain the optimized stencil coefficients. For the other rows the 
regular finite element matrix elements are used. We expect that the fact that no optimized 
coefficients are used inside the PML layer is of little importance, because the accuracy of 
the phase speeds is most relevant for long distance wave propagation, and there is no such 
propagation in the PML layer due to the damping. 

With this method a number of numerical experiments was carried out, varying sev- 
eral parameters. Two choices of medium were considered, the constant and the random 
medium, both on the unit square with PML layers added outside it. A point source at 
(0.3, 0.3) was used a right hand side. We experimented with different multigrid levels from 
2 to 4, using nx x ny= 200 x 200 grid points at the coarse level, and with different sizes 
for the twogrid method, using nx = 200, 400 and 800. Each of the cases was carried out 
with G c = 4 and G c = 6 gridpoints per wavelength at the coarse level. The results are 
displayed in Table ^1 We see that the number of iterations is mostly low (< 20 iterations) 
and increases above this only for the experiments with 3 and 4 levels with a constant 
medium and G = 4 (the most difficult case). 

Hence the new coarse grid operators can be used with PML boundary layers as well, 
providing the same kind of reduction in G c as for the examples without PML. 
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(b) solution with nx = 511, G = 4, a = 2.5e-3 
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Table 5. Convergence results for different methods for a random medium 
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Table 6. Convergence results for different methods for the Marmousi example 
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Table 7. Multigrid convergence for the constant medium with 2,3 and 4 
levels, and different values of G and a 
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Table 8. Two-grid convergence for a constant coefficient medium in 3-D 
using 15 point and 27 point optimized operators at the coarse level. 
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